  function [celltab,xsim]=plot_gamma(alph,bet,mingr,stepgr,maxgr); 
% function [celltab,xsim]=plot_gamma(alph,bet,mingr,stepgr,maxgr); 
% AJ 5/13/04 
% Reburfished version of poltgamma.m 
% See m-file see-bg for description of inputs 
% _________________________________________________________________
meant=alph*bet; 
vart=alph*(bet^2);
sd=sqrt(vart);
if nargin < 3 | isempty(mingr)==1 
    mingr=0.01; 
    maxgr=6; 
    stepgr=0.01; 
else
    if any( [mingr stepgr maxgr] < 0 ) == 1 
        error('Optional Outputs must be positive') 
    end 
end
x=mingr:stepgr:maxgr; 
pdf=gampdf(x,alph,bet); 
figure; 
plot(x,pdf); 
xmm=max(0.01,mingr-2*(stepgr)); 
xlim([xmm maxgr+2*(stepgr)]);
astr=num2str(alph); 
bstr=num2str(bet); 
mstr=num2str(meant); 
sstr=num2str(sd); 
t1=['Gamma(' astr ',' bstr ')   mean : ' mstr ' sd: ' sstr ]; 
title(t1);
% Simulate 
xsim=gamrnd(alph,bet,20000,1); 
per=[0.001 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.999]; 
perc=gaminv(per,alph,bet); 
jj=1; tstr=length(per); 
celltit=cell(tstr,1); 
celltab=cell(tstr,2);
for jj=1:tstr; 
    str_p=sprintf('%1.3f',per(jj)); 
    str_c=sprintf('%3.5f',perc(jj)); 
    celltit(jj)={[str_p,'   ',str_c]};
    celltab(jj,1)={ str_p }; 
    celltab(jj,2)={ str_c };
end 
legend('density',celltit{1},celltit{2},celltit{3},celltit{4},celltit{5},celltit{6},celltit{7},celltit{8},celltit{9})
celladd=[{'1'} {mstr}; {'2'} {sstr}]; 
celltab=[celladd;celltab]; 
format short; 